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INTRODUCTION 


The development of high-speed computers and sophisticated display devices has 
encouraged the development of advanced algorithms for manipulating and displaying 
multidimensional data. The area of computer-aided geometric modeling (CAGM), in par- 
ticular, has advanced significantly during the past few years. In CAGM, computa- 
tional geometry and computer graphics have combined to give mathematical and graphi- 
cal representations of curves, surfaces, and volumes. A wide variety of applications 
may be found in the design of automobiles and aircraft and in the mathematical repre- 
sentation of various physical phenomena, such as geophysical maps and meteorological 
data. 


In CAGM it is often desirable to interpolate three-dimensional surface data con- 
sisting of two independent variables (x and y) and one dependent variable. A num- 
ber of techniques have been developed for surface interpolation, including Coons and 
Bezier patches and tensor products of Bezier curves, cubic splines, and B-splines 
(ref. 1). A difficulty can arise with these methods, especially the spline methods: 
abrupt changes in the dependent variable of the data may induce artificial or exag- 
gerated hills and valleys in the interpolating surface. 

One way to reduce or eliminate the unwanted hills and valleys in an interpolat- 
ing spline surface is to apply tension to the surface. Applying mathematical tension 
to a spline is analogous to grasping the opposite edges of a membrane and stretching 
the membrane to remove wrinkles. 

This paper discusses two algorithms for interpolating surfaces with spline func- 
tions containing tension parameters. Both algorithms are based on the tensor prod- 
ucts of univariate rational splines (ref. 2). The simplest algorithm, which uses a 
single tension parameter for the entire surface, is the birational spline algorithm 
developed by Spath (ref. 2). This algorithm is generalized in this paper to use a 
separate tension parameter for each subinterval along both the x- and y-axes. The 
new algorithm allows for local control of tension on the interpolating surface. Both 
algorithms are illustrated and the results are compared with the results of bicubic 
spline and linear interpolation of terrain elevation data. 


SYMBOLS 

4 by 4 matrix of coefficients 

a i'k£ coefficients of multivariate rational spline in the subregion R^j 

3 (k,& = 1, 2, 3, 4) 

CX^ coefficients in equation for FX ij 

CYj coefficients in equation for FY ij 

coefficients in univariate rational spline on interval i (k = 1, 2, 3, 4) 

dx^ difference between consecutive values of x equal to - x^ 



dyj difference between consecutive values of y equal to y_. +1 _ y . 


F ij 


data value at 


(x., yj > 


FXij 

partial derivative 

of 

fij (x,y) 

with 

respect 

to 

X 

evaluated 

at 

(x i'yj ) 

FY ij 

partial derivative 

of 

fij (x,y) 

with 

respect 

to 

y 

evaluated 

at 

(x^y-j) 

FXYij 

partial derivative 
(x ± ,yj ) 

of 

fij (x,y) 

with 

respect 

to 

X 

and y evaluated at 

fij (x,y) 

rational spline on 

the 

subregion 

R ij 








G i 


g ik (x) 

H j 

hj £ (y) 

m 

n 

P 

Pi 


R ij 


4 by 4 matrix of functions g ik (x) 


functions of x used in rational spline representation on interval i 
(k = 1, 2, 3, 4) 

4 by 4 matrix of functions hj^fy) 

functions of y used in rational spline representation on interval i 
(H = 1, 2, 3, 4) 

number of data points along x-axis 
number of data points along y-axis 
tension factor for entire surface 
tension factor for interval i 
tension factor for interval j 

rectangular subregion in xy-plane defined by x < x < x and 
y . < y < x i+1 


s ij 


4 by 4 matrix of function and derivative values 
r,s,t,u variables used to define rational spline 
x independent variable 

y independent variable 

z dependent variable in univariate rational spline 

^x'^y differences with respect to x and y variables f respectively 
Subscripts: 

i index along x-axis 

j index along y-axis 

k,& general indices 
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Superscripts : 


T matrix transpose 

-1 matrix inverse 

A prime indicates first derivative with respect to the independent variable. A 
double prime indicates second derivative with respect to the independent variable. 


RATIONAL SPLINES 

In this section both the univariate and the bivariate rational splines are de- 
scribed. Because the bivariate rational spline is the tensor product of two univari- 
ate rational splines, the characteristics of the univariate rational spline are 
discussed first. 


Univariate Rational Spline 

Let a given set of n data points be represented by (x^,z^), where 
i - 1,2,...,n and x 1 < x 2 < *** < x n « The values of the independent variable x.^ 
need not be equally spaced. The rational spline on interval i (i = 1,2,...,n-1) is 
defined in references 2 and 3 to be 


4 

z = ^ c., g., (x) (x. < x < x ) (1 

*“■ ik ^fk l i+1 

k=1 

where c^ are unknown coefficients, p^ is the tension parameter for interval i, 
and 

3 

gi1 (x) = u g i3 (x) = " + y 

l 


g i2 (x) = t <?i4 (x) = p.u +1 

i 

where 


u 


X i+1 


- x 


dx. 

l 


t 


x 


X. 

1 


1 


u 


dx. 

l 


3 



dx i = 


\L+1 


- x- 


Equation (1 ) is defined for all values of the independent variable x in the 
data range if the tension parameter Pi is restricted to p. > -1. if p . i s set 
to zero, equation (1 ) reduces to a cubic spline function. As p. increases from 
zero, the cubic terms decrease in magnitude and the function tends to the equation of 
the line joining the knots at x^ and x^ +1 . Because a distinct, independent ten- 
sion parameter is associated with each interval, the behavior of the function in each 
interval may be locally controlled. 


Evaluation of equation (1) for each subinterval requires knowledge of the four 
coefficients c^^ , c ^2' c i3' anc ^ c i4* T ^ us f for n data points (equivalently, 

for n - 1 subintervals), 4n - 4 coefficients must be determined. Spath (ref. 2) 
reduces the magnitude of this problem by writing the coefficients in terms of the 
values of the function and first derivative at the knots. End conditions are applied 
to the first derivative, and equations ensuring the continuity of the second deriv- 
ative at the interior knots are derived. This derivation yields a system of n - 2 
equations for the ^n - 2 unknown interior first derivatives. Frost and Kinzel 
(ref. 3) extend Spath' s approach by allowing for three different end conditions and 
by developing an iterative method for determining the tension parameters. Tension 
parameters are found so that the rational spline deviates from the line joining knots 
by no more than a prescribed value. 

Another approach to determining a good fit to the data was taken by Schiess and 
Kerr (ref. 4) in deriving a least— squares rational spline approximation. The ra- 
tional spline is reformulated in terms of the unknown spline function and its second 
derivative at the knots. The conditions imposed at the interior knots are that the 
first derivatives be continuous. This leads to a constrained least— squares problem 
in the 2n values of the unknown function and its second derivative at the knots. 


Bivariate Rational Spline 

Let a given set of m by n data points in three dimensions be represented by 
( X i/Yj/Fij), where i = 1,2,... ,n and j = 1,2,...,m. It is assumed that the inde- 
pendent variables are ordered (x 1 < x 2 < ••• < x n and y., < y 2 < • • • < y ) and form 
a rectangular grid, but are not necessarily equally spaced. An interpolating surface 
through the given points is desired. 


(i 


The bivariate rational spline on the subregion R. . defined by x. < x < x 
1,2,...,n-1) and < y < Yj +1 (j = 1#2,.,.,m-1) is defined in reference 2 + by 


fij (x/Y) 


4 4 

-L E 

k=1 £=1 


a ijk& 


g ik (x> V y) 


( 2 ) 
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are 


where g^(x) and p^ are the same as for the univariate spline, a ijk£ 
unknown coefficients, qj is the tension parameter for interval j, and 


h jl( y)-s h j 3 W - q. r + 1 


h j2 (y) = r 


V (y) = 5Ts~rT 


where 


Vi - y 

S = —*~z 

dy. 


r = 


y _iZi 


= 1 - s 


dy j " y i + 1 - y j 


As defined by equation (2), the bivariate rational spline on each subregion R. . 
is a function of 2 tension parameters and 16 coefficients. The coefficients are 
determined so that the rational spline and its first and second derivatives are con- 
tinuous over the entire region. The m + n - 2 tension parameters may be adjusted 
individually; each parameter affects the behavior of the rational spline in a strip 
parallel to either the x-axis (for q j ) or the y-axis (for p^). In this paper, this 
form of rational spline is called the J multiple-parameter rational spline. 

Since 16 coefficients are needed on each subregion, a total of 16(m - 1)(n - 1) 
coefficients must be determined to define the entire rational spline. For example, 
for a 30 by 30 grid (m = n = 30), 13 456 coefficients are needed. Spath (ref. 2) 
reduces the actual number of unknown quantities by writing the coefficients as linear 
combinations of the values of the function and its derivatives at the grid points. 

Let and FY^j be the first derivatives of f^j(x,y) with respect to x 

and y, respectively, and FXY^j the cross derivative, all evaluated at the point 
(Xi, yj ). For the subregion j , define the following 4 by 4 matrices: 


S ij 


' F ij 

FYij 

F i ( j+1 ) 

FY i(j+1) 

Fx ij 

FXYij 

FX i( j+1 ) 

FXY i(j+1> 

F (i+1 ) j 

FY (i+1)j 

F (i+1)(j+1) 

FY ( i+1 ) (j+1 ) 

_ FX (i+1 ) j 

FXY (i+1)j 

FX (i+1 ) (j+1 ) 

FXY (i+1 ) ( j+1 ) 
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a i j 1 1 

a i j 1 2 

a i j 1 3 

a ij 1 4 

a i j 21 

a i j 22 

a i j 23 

a i j 24 

a i j 31 

a i j 32 

a i j 33 

a i j 34 

_ a ij41 

a i j 42 

a i j 43 

a i j 44 


g i1 (x i ] 

g i 2 ^ x i ^ 

g 13 (x i> 

g 14< x i> 

g! . (x. ) 

il i 

g l 2 <X i ) 

g i3 <X i> 

g i4 (x i> 

g i1 ^ x i+1 

> 9i2 (x i+l> 

9i3 (x i+1> 

g i4 < x i + i) 

_ g ii {x i+ i 

) g: 2 (x i+1 ) 

g i3 (X i + 1> 

g i4 (x i+l>_ 

“ 1 

0 1 

0 


1 

dx. 

l 

1 3 + p 

dx. dx, 

l l 

^ 0 


0 

1 0 

1 


1 

1 0 

3 + p. 

l 


dx. 

dx. 

dx. 


_ l 

l 

i _ 


“ h jl (y j ) 

hj 2 (y j ) 

h J3 (yj) 

h j4 (yj)- 

W 

h j2 ( V 

W 

w 

h ji (y j+i 

> h j2 ( yj+i> 

h :3 (y j+i > 

h j 4 (y j+1 ) 


h j2 <y j+1 ! 

V y j+ 1> 

h j4 (y j+ ,L 

1 0 1 

i i 3 + 

dy. dy. dy. 

0 

0 


0 

1 0 

1 


1 

J- 

3 + q . 
3 


L dy j dy j 

dy . J 
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Then in matrix notation 


S- • = G. A. ,H T (3) 

i: i 13 3 


Equation (3) can be verified by differentiating equation (2), as appropriate, and 
evaluating the results at the corners of the subregion. Note that and Hj 

depend on the grid spacing and tension parameters but not on the function values. 

Since the matrices G^ and Hj are nonsingular, equation (3) can be solved for 
the matrix of coefficients: 


A ij 




-1 


(4) 


Therefore, for any subregion the 16 coefficients can be determined from the values of 
the function, its first derivatives with respect to x and y, and its cross deriv- 
atives at the four corners of the subregion. Therefore a total of 4mn function and 
derivative values are needed to calculate the coefficients. In addition to the given 
function values, the algorithms to be described require values for the derivatives 
FX^j on the boundaries x = x 1 and x = x n , for the derivatives FYqj on the 
boundaries y = y^ and y - y m , and for the cross derivatives at the four corners 
of the region. This reduces the problem to one of solving for 3mn - 2(m + n) - 4 
derivative values. The magnitude of this problem is further reduced by decomposing 
it into four smaller subproblems, each of which consists of solving tridiagonal 
systems of equations. 


RATIONAL SPLINE ALGORITHMS 

In this section algorithms for finding the surface-interpolating rational spline 
in terms of the function and its derivatives are presented. The first algorithm pre- 
sented is for the general case of n — 1 values of the tension parameters p^ 

(i = 1,2,...,n-1) and m - 1 values of the tension parameters qj (j = 1 ,2, . . . ,m-1 ) . 
This is a new algorithm not presented elsewhere; the algorithm is J derived in the 
appendix. The second algorithm is for the special case of a single tension param- 
eter P for the entire surface. This algorithm was originally developed by Spath 
(ref. 2). 


Multiple-Tension-Parameter Algorithm 

Let the data (x^,yj ,F^j ) and tension parameters p^ and qj (for 

i = 1,2,... ,n and j = 1,2,...,m) be given. In addition, 2m + 2n + 4 boundary 

conditions must be given: the derivatives with respect to x along the boundaries 

x = x^ and x n ( FX ij and FX n j ' 3 = 1/2,.. ,,m), the derivatives with respect to 

y along the boundaries y = y 1 and y m (FY^ and FY^ m , a = 1/2,..., n), and the 

cross derivatives at the corners (FXY 1 1 , FXY. , FXY . , and FXY ). 

11' 1m' nl ' nm 
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The algorithm for finding the remaining 3mn - 2(m + n) — 4 derivatives pro 
ceeds in four stages: 

1 . Find the derivative with respect to x at each interior point; for each 
value of j (j = 1 f 2,...,m), solve for n - 2 values of FX^j in the equations 


CX i-1 FX ( i— 1 )i + [ ( 2 + P, - )CX. + (2 + p. )CX. ]FX. . + CX. 

' ' ■> l—i i i i 


FX.. , X- 
(i+l )d 


3 + P i-1 
dx i-1 


CX, , A F,, + 


3 + p, 


-v. ■ , CX. A F. . 

1-1 x (i— 1 )] dx. l x ij 


(i = 2,3,.. . ,n-1 ) 


where 


P i + 3 ^i + 3 

[ (2 + p ^) 3 - 1 ] dx ^ 


CX, = 


A x F ij F (i+1)j F ij 


2. Find the derivative with respect to y at each interior point; for each 
value of i (i = 1,2,...,n), solve for m - 2 values of FY^ j in the equations 

CY j-1 FY i(j-1) + [ <2 + q j-1 ,Cr j-1 + 12 + qjIC^lFY.. + CY. FY 1( . +1) 


3 + q j-1 3 + q 

— ; CY. A F. . + 2. CY. A F 

dYj.-, D-1 Y i(D-1) dy j y r ij 


( j = 2,3,... ,m-1 ) 


where 



3. Find the cross derivatives along the boundaries y = y., and y m ; for j = 1 
and m, solve for n - 2 values of FXY^ in the equations 


CX 


i -1 


FXY 


(i -1 )j 


[(2 


P. Jcx. 
1-1 1-1 


+ (2 + p. )CX. I fXY. . + CX. 

l l J 13 1 


FXY 


( i +1 )j 


3 + p 


i -1 


dX i -1 


CX. A FY,. . . 
1-1 x (l- 1)3 


3 + P i 

+ — =- CX. A FY. . 

dx. 1 x 13 


(i — 2,3,.../ n— 1 ) 


(7) 


where 


A FY. • 
X 13 


" FY (i+ 1 )j “ FY ij 


4. Find the cross derivatives at the interior points; for each value of i 
(i = 1 , 2 , . . . ,n ) , solve for ra - 2 values of FXY^ in the equations 

C*j-1 FXY i(j-1) + [< 2 + - < 2 + + CY j FXY i(j + 1 ) 


3 + q i-1 3 + q j 

S_L cy. , A FX. , . „ , + — ; — CY. A FX. . 

dy^ 3-1 y 1 ( 3 - 1 ) d Yj i y 13 


( j = 2,3, .. . ,m-1 ) (8) 


where 


V x ij 


= FX 


i ( j+1 ) 


- FX 
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Equations (5) to (8) have a similar form; in fact the coefficients of equa- 
tions (5) and (6) are identical to those of equations (7) and (8), respectively. The 
similarity of the equations can be exploited to reduce computational requirements 
since only two distinct matrices must be inverted. 

The final step of the procedure is the application of equation (4) to the deriv- 
atives in order to calculate the coefficients needed for equation (2). Depending on 
the availability of computer time and memory, the user may prefer to store the deriv- 
atives and to calculate coefficients only as they are needed. 


9 



Single-Tension-Parameter Algorithm 


The algorithm requiring only one tension parameter value P uses the same data 
and derivative values that are required for the multiple-tension-parameter algorithm. 
However, the equations for this algorithm can be obtained from the equations for the 
preceding algorithm by setting all the tension parameters equal to P (p. = p f or 
i 1,2,...,n-1 and qj - P for j = 1 , 2 , . . . ,m-1 ) . Using a single tension param- 
eter simplifies the coefficients of equations (5) to (8) considerably. 

The four stages of the algorithm are as follows: 

1. Find the derivative with respect to x at each interior point; for each 
value of j (j = 1,2,...,m), solve for n - 2 values of FX^ in the equations 


— "u-Dj + (2 + p »(ar 


dx i-i 


i-1 


+ V x - • + T" — FX,. 
dx i J 1 J dx i ( 1 + 1)3 


fh F. 


A F. 


(3 + P)f- X + -2-iI 


dx Li 


dx. 


(i = 2,3, .. . ,n-1 ) 


(9) 


2. Find the derivative with respect to y at each interior point; for each 
value of i (i = 1 , 2, . . . ,n) , solve for m - 2 values of FY — in the equations 


FY. 


dy. T Kj-1) 


+ (2 + P )(~ + W. . + FY. 


dy. 


dyl ij dy i ( j+1 ) 


/A F 

{3 + P ) f . y 1(3 ~ 1} + y 


dy 


j-i 



( j = 2,3,... ,m-1 ) 


( 10 ) 


3. Find the cross derivatives along the boundaries y = y anc j y . f or j = 1 


and m, solve for n — 2 values of FXY^j in the equations 


dx FXY fi-1H + ^ 2 + p ) f + X~“V XY • • + ^ — FXY , . , 

1-1 U 1)3 \ d i-1 dx i/ 13 dx i <i+1 )j 


M FY.. 
" _U 

2 


(3 + p)l22LlSlzlll + ,Ji FY d3 


dx 


i-1 


dx. 


(i = 2,3, .. . ,n-1 ) 


(11 ) 


10 



4. Find the cross derivatives at the interior points; for each value of i 
(i = 1,2,...,n), solve for m - 2 values of FXY — in the equations 


FXY ... . , + (2 + P ) ( ^ + AFXY . . + — — FXY . . . . . 

dy_._ 1 i( 3-D \ y j-1 dy j/ 13 dy j l( 3+D 


/A FX.,. „ . A FX. .' 
(3 + P)f Y ■ 2 3 ~ - + - JL ~ P 

dy j-i a h 


( j = 2 , 3 , ... ,m-1 ) 


( 12 ) 


Equations (9) to (12) are identical to equations (8.25) to (8.28), respectively, 
of reference 2. All the equations have the same general form and the coefficients of 
equations (9) and (10) are identical to the coefficients of equations (11) and (12), 
respectively. 


Selection of Tension Parameters 

Although an automated procedure for adjusting tension parameters is available 
for univariate rational splines (refs. 3 and 4), no such procedure has been devised 
for bivariate rational splines. Instead, it is recommended that the user visually 
examine three-dimensional or contour plots of the data and the rational spline inter- 
polating surface and then use engineering judgment to select tension-parameter 
values. This approach is outlined here. 

First the user should obtain a plot of the original data. This plot will show 
both general trends and local anomalies of the data. Second the user should compute 
the interpolating cubic spline surface by calculating the interpolating rational 
spline with all tension parameters set to zero. Comparison of the data and the cubic 
spline plots will indicate regions where the cubic spline exhibits undesirable or ex- 
aggerated hills and valleys. Using small to moderate tension-parameter values (say, 
from 1 to 10) for those regions, the user can recalculate an interpolating rational 
spline. In this way, after a few trials with different tension values, a rational 
spline surface can be obtained which the user considers representative of the data. 


EXAMPLE 


For this example a set of data was chosen to show the flexibility and capabili- 
ties of the interpolating rational spline. The data consist of terrain elevation 
measurements on a rectangular area of size 2600 by 2000 feet. The measurements were 
taken at 100-foot intervals in both the x- and the y-direction. Elevations were 
measured to an accuracy of 0.1 foot. 

Figure 1 shows a surface plot of the terrain elevation data. Note the cliff 
along the line x = 1800 feet. Because of this cliff, extra measurements were taken 
at the base of the cliff along the line x = 1795 feet. Including this set of mea- 
surements, the size of the data grid is 28 by 21. The elevations range from 7.9 feet 
(at x = 2300 ft, y = 100 ft) to 42.6 feet (on the cliff at x = 1800 ft, 
y = 800 ft). The maximum rate of change in elevation also occurs at the point of 
maximum elevation, where the elevation increases 11.2 feet in a distance of 5 feet 
(from x = 1795 ft to x = 1800 ft). 
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It is desired to produce a finer grid by interpolating the data every 25 feet 
in both the x- and the y-direction. As a basis of comparison with the rational 
splines, a bilinear interpolation was first performed. The bilinear interpolation 
consists of a linear Lagrange interpolation in each of the x and y variables 
(ref. 5). Figure 2 shows a surface plot of the bilinear interpolation of the terrain 
data. 


Computer software used to determine the interpolating bivariate splines in this 
study is based on the software published by Spath (ref. 2). For the single-parameter 
rational spline the FORTRAN code given in reference 2 was used. This code solves 
equations (9) to (12). For the multiple-parameter bivariate rational spline the 
Spath software was modified to solve the multiple-tension-parameter equations 
(eqs. (5) to (8)). Neither the original Spath software nor the modified software has 
been optimized to take advantage of the similarity in the respective sets of equa- 
tions. Rather, each of the equations is re-solved as many times as are dictated by 
the algorithms. Timing comparisons were made on a Control Data Corporation CYBER 175 
computer at Langley Research Center. For the example given here the required central 
processor times were 4.7 seconds for the bilinear interpolation, 9.2 seconds for the 
single-parameter rational spline, and 9.4 seconds for the multiple-parameter rational 
spline. Both rational spline algorithms therefore require approximately the same 
amount of central processor time. Times for both rational spline algorithms can be 
reduced considerably by taking advantage of the similarity of the equations solved. 

Derivative information along the surface boundaries is needed for both rational 
spline algorithms. A simple way to estimate these derivatives, which is also the 
method used in the example, is to calculate one-sided finite differences. These dif— 
f®^®^ces are calculated from data points on the boundaries and immediately interior 
to the boundaries. 

With the tension parameters set to zero in either bivariate rational spline 
function, a bicubic spline function results. Figure 3 illustrates the interpolating 
bicubic spline surface for the terrain data interpolated at 25-foot increments in 
each direction. For univariate data having a drastic change in slope, the cubic 
spline typically deviates extensively from the desired trend between data points 
(ref. 3). Figure 3 shows that the same phenomenon occurs with the bicubic spline. 

The interpolating bicubic spline has exaggerated the height of the cliff between data 
points to an elevation of approximately 75 feet (at y = 925 ft). At the same time 
another oscillation in the surface has created a valley and hill in front of the 

In some areas the bicubic spline surface (fig. 3) appears to be slightly more 
rounded with low hills where the bilinear surface (fig. 2) is essentially flat. 

The single-parameter rational spline was applied first in order to reduce the 
artificial hills and valleys created by the bicubic spline. Recall that this type of 
rational spline is analogous to grasping the surface at all the edges and pulling to 
stretch the surface. Figure 4 shows the result of applying a tension of 40. Most of 
the exaggerated cliff height and the artificial valley and hill in front of the cliff 
have been eliminated. The surface in this figure closely resembles the original 
surface in figure 1 . 

The advantage of the multiple-parameter rational spline is that it allows local 
control of surface undulations through selective assignment of nonzero tension param- 
eters. In this example it is especially desirable to eliminate surface oscillations 
in the neighborhood of the cliff. This can be accomplished by applying nonzero 
tension in the neighboring intervals parallel to the cliff. A tension value of 40 
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2600 



Figure 4.- Single-parameter rational spline interpolation of terrain 
elevation data with tension of 40. 


2600 



Figure 5.- Multiple-parameter rational spline interpolation of terrain 
elevation data with tension of 40 for 1700 ft < x < 1900 ft; i.e., 

p 18 = P 1 9 = p 20 = 40 * 


2600 



2000 


Figure 6.- Multiple-parameter rational spline interpolation of terrain 
elevation data with tension of 100 for 1700 ft < x < 1900 ft; i.e., 

P 1 8 = P 1 9 = p 20 = 100 ‘ 
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(which was used for the single-parameter rational spline) was assigned to param- 
eters P-|g' Pi g, and P20* T * ie nonzero tension parameters are in the range 
1700 ft < x < 1900 ft, which is centered about the cliff ridge. All other tension 
parameters are zero. The result of this rational spline interpolation is shown in 
figure 5. This figure shows that the cliff elevation has been reduced considerably 
and the valley and hill in front of the cliff eliminated. The remaining areas of 
this surface are the same as the bicubic spline surface. 

As a further test of the multiple-parameter rational spline, the values of the 
same three nonzero tension parameters were increased to 100. Figure 6 shows that the 
larger tension reduces the cliff height even further. Note, however, that most of 
the excessive height was reduced by increasing the tension from zero to 40. The 
result of the tension value of 100 is to produce an interpolating surface that is 
essentially linear along the x-direction near the cliff and cubic in x and y 
otherwise . 


CONCLUDING REMARKS 

Two algorithms for surface interpolation with bivariate rational spline func- 
tions on a rectangular grid have been presented. The rational spline functions com- 
bine the advantages of a cubic function having continuous first and second deriva- 
tives throughout the interpolatory region with the advantage of a function having 
variable tension. Adjustment of one tension parameter in the single-parameter 
rational spline allows the user to reduce unwanted oscillations to any desired extent 
across the entire surface. The multiple-parameter rational spline provides adjust- 
able parameters for each rectangular subregion and thus gives the user control over 
local behavior of the interpolatory surface. A new algorithm for finding the coeffi- 
cients of the multiple-parameter rational spline has been derived and presented. 

The terrain elevation example presented illustrates the reduction in undesirable 
oscillations that is possible with a bivariate rational spline. The example demon- 
strates that the single-parameter rational spline can be adjusted to behave like any 
interpolatory surface ranging from a bilinear to a bicubic spline surface. The exam- 
ple also illustrates that the multiple-parameter rational spline provides the cap- 
ability of controlling local oscillations in the surface caused by abrupt changes in 
trends in the data. 

There are two major disadvantages to using the bivariate rational spline. 

First, the algorithms for finding the rational spline coefficients and interpolating 
points may require significantly more central processor time than does a bilinear 
interpolation. However optimizing the computer code for these algorithms can reduce 
the necessary central processor time. Second, several computer runs with different 
tension values may be necessary to find an interpolatory surface satisfactory to the 
user. The single-parameter rational spline in the terrain example required four 
computer runs to attain the results presented. These results, however, immediately 

led to the first results of the multiple-parameter rational spline interpolation. 

The second disadvantage can be overcome with experience in applying the rational 
splines. An alternative approach would be to apply an expert system to selection of 
the tension parameters. 


NASA Langley Research Center 
Hampton, VA 23665-5225 
October 28, 1985 
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APPENDIX 


DERIVATION OF CONTINUITY EQUATIONS 

In this appendix equations (5) to (8), which are solved for the unknown first 
derivatives and cross derivatives at interior grid points, are derived. All four 
sets of equations have comparable derivations. 


Derivation of Equation (5) 

Using the definition of the rational spline in equation (2), consider the 
interval x.^ < x < x i+1 , fix j, and evaluate f ij(x,y) at y . . Then 
h j i ( Yj } = h j3 (yj) = 1, h j2 (yj) = h j4 ( Yj ) = 0, and 


f ij(x, yj ) 


4 

L 

k=1 


g (x) (a. 
lk ljkl 


a. ., _ ) 
i]k3 


(A1 ) 


For simplicity, define b^j^ = a ijk1 + a ijk3* T h en from equation (A1 ) and the defi- 
nitions of ^i^x), 


f ij(Xi,yj) = b ±j1 + b ij3 = F 




f ij (x i+1' y j) b i j 2 + b ij4 " F (i+1 )j J 

where F — and F (i + -j)j are known function values at the grid points, 
The first derivatives of g ik (x) are 


(A2 ) 


5 il <x) 


dx. 

i 


g' 2 (x) 


dx. 

l 


g i3 <x) ‘ 


g: 4 (x) - 


-3u 3 (p^t + 1) - u 3 p, 

dx . (p . t + 1 ) 3 
1 1 


3t^(p^u + 1) + t^p^ 


dx. (p.u + 1 ) 

l l 




J 


(A3) 
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Differentiating equation (A1 ) with respect to x, applying equations (A3), and evalu- 
ating at x^ and x^ + .j yields 


3f. . 
3-3 


(x. ,y. ) = - 


/V \ . f 1 • 

dX 1 J 


b b 3 + p. 

in 1 ij2 i . 

— - — — + — — - b. ._ = FX . . 

dx. dx. dx. ij3 lj 




3f 


3-3 


(x. , wY.) = " ' 


3x ' i+1 ,jr j 


b. . b. 3 + p. 

ill . 3.J2 1 . = py 

dx_. dx. dx. ij4 (i+1)j 


(A4) 


where FX^ - and FX (i+i)j are the unknown first derivatives with respect to x at 
interior grid points. 

The second derivatives of g i)c (x) are 


(x) = gv 2 (x) = 0 


g i3 (x) " 


g^ 4 (x) = 


6u(p^t + I) 2 + 6u 2 p^(p^t + 1) + 2u 3 p 2 

dx 2 (p. t + 1 ) 3 
l i 


6t(p.u + I) 2 + 6t 2 p.(p.u + 1) + 2t 3 p 2 

l l l l 

dx?(p.u + 1 ) 3 

l l 




J 


(A5) 


Differentiating equation (A1 ) twice with respect to x, applying equations (A5), and 
evaluating the results at x^ and x^ +1 leads to 


3 2 f . . 


3x 


£•<*.. y.> - 


2 i j 


2p . + 6p . + 6 

l l 

dx 2 

l 


b. .. 
3-3 3 




3 2 f . . 
11 


3x 




2 i+1 ' 


2p. + 6p. +6 

1 1 

dx 2 

i 


b ij4 


J 


(A6) 
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Solve equation (A2) for and substitute into equations (A4), 


rearrange the results to obtain 


and 


-(2 + Pl )b lj3 - b ij4 - d Xi FX ij + Py - P (i+1)j 


(A7 ) 


b ij3 + < 2 + Pi> b ij4 * d *i ra (i+1)j + Fij - F( 1+1)j 


( A8 ) 


Define - F (i+i ) j - F ij and solve equation (A7) for b ± j 4 to get 


b ij4 - -(2 + P i )b. j3 - d X;L FX ±j + A x F ± . 


( A9 ) 


Substituting equation (A9) into equation (A8) and solving for b^j^ yields 


13 + h’Vn ~ d *i FX (1+1 )i - + P^x. PX 14 


i j 3 


(2 + p . ) - 1 

l 


(A10) 


Substituting equation (A10) into equation (A9) gives 


-(3 + p )A F + (2 + p )dx FX. . . . + dx. FX. . 

b. . . = — 2 - JU l l (i +1 n 3 . l] 

ID 4 2 

(2 + p. ) - 1 

1 


(All ) 


The continuity condition to be imposed at x. is that the second derivative is 
continuous: 


f (i-1 )j <X i’ y j * " f h ( VV 


or from equations (A6) f 


p i-i + 3 p i-i + 3 


(i-1 )j4 


p. + 3p. + 3 
dx 2 " 

l 


b. 

ID 3 
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Substituting equation (A10) and equation (All) with i replaced by i - 1 and rear- 
ranging yields 


.2 (i-1)D 


[ ( 2 + Pi.-, > " 1 ] dx i _ 1 


(p?., + ^i.T + 3 ) (2 + P^) 
[(2 + P^,) 2 - 


( p i + 3p ± + 3) (2 + p i ) 

[(2 + p.) 2 - l]dx. J 1D [(2 + p.) 2 - ijdx. 


p . + 3p . + 3 
1 1 

>FX. . + FX,. ... 

2.1, (l+l )] 


(p , + 3p. , + 3) (3 + p. , ) 
v 1-1 1-1 J 1-1 


( p i + 3p i + 3) (3 + p ± ) 


[ (2 + Pi_, ) - l]dXi_, 


A F + = — A F. . 

2 .1,2 x (i -1 ) j F 0 ± „ \2 _ ,i x 1 j 


( A 1 2 ) 


[(2 + P . r - i]dxf 

L ^1 J 1 


Equation (A12) is identical to equation (5), which holds for j = 1,2,...,m and 
1 = 2,3,4.. ,n-1 . 


Derivation of Equation (6) 

The derivation of equation (6) is very similar to the derivation of equa- 
tion (5). For this derivation, fix i and consider the interval yj < y < yj+-j* 
Proceed through the same steps as before, but consider the first and second deriva- 
tives with respect to y of f^j(x^,y) evaluated at yj and Yj + i* Equations (5) 
and (6) have a similar form because the derivations are entirely analogous. 


Derivation of Equation ( 7 ) 

The derivations of equations (7) and (8) are very similar to the derivations of 
equations (5) and (6). The major difference is that the conditions imposed in this 
derivation and the next are the continuity of the mixed third-order derivatives at 
the grid points. 
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First consider the first derivative of equation ( 2 ) with respect to y evalu- 
ated at y = Yj, which will remain fixed. This derivative can be written 


^ = h h 9ik<x) w 


4 

E 

k=1 


a ijk1 a ijk2 3 + q j \ 

^ dy_. dy_. dy_. a ijkl) g ik^ X ^ 


4 

E 

k=1 


d. g (x) 
13k lk 


(A 13 ) 


where d£j k has been defined for convenience. Differentiate equation (A 13 ) twice 
with respect to x and evaluate at x^ and x i+1 to get 


3 Vvy 

9 2 x 9y 


d. 

13 3 


2p . + 6p , + 6 

1 

dx 2 

1 


8 f ii (x i + 1- y j ) 
9 2 x 9y 


d ij 4 


2p. + 6p. +6 

1 1 

dx 2 

1 




> 




<A 14 ) 


The continuity condition to be fulfilled is that the third-order mixed deriva- 
tives are continuous at x^: 


3 Vih ( w 

9 2 x 9y 


f . . (x. ,y . ) 
9 2 x 9y 


or from equations (A 14 ), 


P i-1 + 3p i-1 + 3 
(i-1 )j 4 2 

dx i-1 


d ij 3 


p . + 3 p . + 3 
1 1 

dx 2 

1 


(A 15 ) 
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Equation (A15) may be used to solve for FXY^j (i = 1,2,...,m). In order to do so, 
the derivative with respect to x of equation (A13) and equation (A13) itself must 
be evaluated at x^ and x^ +1 . Evaluating equation (A13) at x^ and x ^ +1 gives 


3f . . (x. ,y . ) 
il 3- 3 

3y 


3f . . (x ,y . ) 

13 1+1 3 


3y 


d + d. = FY. . 
13 I 13 3 13 


d. . _ + d. . . = FY.. . . . 
132 134 (l+l >3 


(A16) 


J 


Solving equations (A16) for d.^ and d^j 2 in terms of d^-j, d ij 4 » FY i j ' and 
FY (i + 1 )j Yields 


d. . . = FY- • - d. • 
13 1 13 133 


d ij 2 FY (i+1)j d ij4 




(A17) 


Since equation ( 6 ) was previously solved for FY ij (i = 1,2,...,m; j = 1,2,...,n), 
these quantities are now known. Differentiating equation (A13) with respect to x 
and evaluating at x i and x i+1 gives 

'N 


> (A18) 


J 


3 f . . (x. ,y . ) 

13 1 3 _ 1 


3x 3y 


[-d. + d. - (3 + p. )d. - FXY. . 

dx i L 13 1 132 1 i]3 J 13 


13 


3 f . . (x ,y. ) 

13 1+1 3 _ 1 


3x 3y 


[-d. + d. + (3 + p. )d. . ] = FXY,. 

dx^ L ljl 132 1 134 j (i 


j4 J *“‘(i+1)j 


Substitute in equations (A18) for and d ij2 f rom equations (A17) and rear- 

range to get 


-(2 + Pi )d ij3 - d . j4 = d Xi FXY i;j + FY i;j - PY (i+1)j 

> 

d i j 3 + < 2 + Pi^ d ij4 = dx i FXY (i + 1 )j + FY ij " FY (i+ 1 )j 


(A19) 
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Using the notation A x FY^. = FY,^ + ^ v. - FY. ■ , equations (A19) can be solved for d., • o 
and d ij4 to yield J J J J J 


(2 + p. )dx. FXY. . + dx. FXY, . - (3 + p )A FY 

d _ i l il l ( 1 +D 3 i x ij 

ij3 .2 


A 


1 - (2 + p. ) 

l 


(A20) 


M. j 4 


-dx. FXY. . - (2 + p. )dx. FXY, . , „ , . + (3 + p, )A FY 
i 13 i i (i+ 1)3 1 x 


13 


1 - (2 + p. )' 

1 


J 


Use equations (A20) to substitute for d^ and d (i-1) j 4 in equation (A15); then 
combine and rearrange terms so that the unknown cross derivatives are on the left 
side of the equation. The final result is 


h-i + 3p i-i + 3 


<A-i + 3p i-1 + 3 1< 2 + Pi., I 


— FXY + ^ 

dx-.jl - (2 +P._ 1 ) ] U ' ~' 13 (_ dx w [l - (2 + Pl-1 ) 


( p i + 3p ± + 3) (2 + pj 


dx , [ 1 - (2 + P . r ] 

i L 1 J 


p . + 3p . + 3 

n , FXY. . + - - — FXY. 

2 ^ ( 13 


dx. [l - (2 + p. ) 2 ] (l+1 )j 

l 1 J 


l p i-i - 3p i-, + 3 ) (3 +p i-,> 

*<■_,[’ - (2 + p.^) 2 ] 


A FY, . ... 
x (1-1)3 


(p? + 3p i + 3) (3 + P ± ) 

dx 2 [l - (2 + p. ) 2 
1 L 1 


-A FY. . (A21 ) 

x 13 


Equation (A21 ) is solved for FXY^ for i = 2,3,...,n-1 and for j = 1 and m; 
this gives the cross derivatives along the boundaries y = y 1 and y m . 

Derivation of Equation (8) 


Equation (8) is derived in a manner similar to the way equation (7) is derived, 
However for this derivation the continuity condition 


9 3 f . 


3 

9 f. . 


i ( j — 1 ) , x ij , 

— i (x. ,y. ) = (x, , y ) 


9x3 y 


3 x 3 2 y 1 ” ;I 
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is imposed. This is accomplished by fixing i and using an expression analogous to 
equation (A13) for the first derivative with respect to x. The derivation then 
proceeds comparably to the derivation of equation (7). 
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